/***************************************************************
But still has a problem where starting with the same seed clone does
not give the same results -- where it does for agarose.
                             cb_assemble
previously called  zcalc.c 
Contains routines to calculate cbmaps - used in cb_compute_order.c and cb_known_order.c

bug fix: W. Nelson
Problem was that the CB bands for a clone were stored in an array of size MAX_BANDS.
Since some of the CBs don't actually match to the clone, there could be more of them 
than the number of bands in the clone. If the clone had alot of bands to begin with, 
the total of CBs could exceed MAX_BANDS. Solution was to use twice MAX_BANDS for 
those arrays. 

CARI 10.feb.6 = add logic so that HICF uses percentage of gap.
                and agarose uses original values
***************************************************************/
#include <stdio.h>
#include <malloc.h>
#include <math.h>
#include "clam.h"

int QTRACE=0;
extern int tmptest;

/* these are for agarose, so algo works like use to */
#define AG_GAP 3
#define HICF_GAP 8
#define AG_HANG_EXTENT 10
#define HICF_HANG_EXTENT 15

#define LEFT 0
#define RIGHT 1
#define CENTER 2
#define DUD 3
#define ALIGN 4
#define PQ_BURIED 50002

#define Q_RULE(a) ((int) ((float) a - ((float) a * Pz.p_align)+0.5))
#define END_RULE(a) (a * 0.3)
          /* number of extra set bands to check against fp */
#define WIN_FUZZ_RULE(a) (floor(a * 0.2 + 0.5))
#define WIN_BURY_FUZZ(a) (a * Pz.burFlt)
#define GAP(a) (floor(a * Pz.gap + 0.5))

#define BAD(msg) {printf("%s\n",msg); Zprepare_for_display(); Zclamdump();  return;}
#define BAD1(msg) {printf("%s\n",msg); Zprepare_for_display(); Zclamdump(); return 0;}

extern int  Zget_tol();
extern void loadCz(), Zprepare_for_display(), Zdump_matrix(), Zdump_all(), Zclamdump();
extern int Zcreate_node();

int Zcnt_extend; /* used in shuf to know when to shuf */
int Zbcm(), compare_to_Zset();
int make_this_Zset(), find_Zgroup(), get_Zgroup(), get_Zset(), remove_any_bad_bands();
void init_Zset(), Zprepare_for_display();
void roll_left(), roll_right(), add_to_hang(), extend_Zgroup(), maybe_add_Zgroup();
void add_extra_ZZ(), add_clone_Zset(), add_nodes_to_pq();
ZPQZ *insertZpq(), *drainZpq();
void cbmap_sort();

#define TBAND 3*NBANDS

struct tmpZ {        /* this is externed in zorder.c, so any changes...!! */
   int bit[TBAND], who[TBAND], extra[TBAND];
   int score, match, nomatch, n_extra;
   int leftb, rightb, left_no, right_no;
   int left_sp, right_sp, left_ok, right_ok;
   int fstg, lastg;
   int ok, high, max;
} Tz; 

static void updateZZ(), msg_clone();
static int  get_Zband(), split_Zgroup();

static struct winy {
   int avg, g, who;
} *win=NULL;

static int n_win=0, C, First_rightb, remove_bad=0;

static int Zrule_avg=0;  /* 0 = 1= average ,medium 2=midpt */
static int Zrule_shave=1;  /* do 2xtol*/
static int Zrule_no_1_grp=0; /* no groups of size one - doesn't help */
static int Zrule_remove=1; /* remove cb if o's > +'s */
static int Zrule_pq=1;  /* 0 - exponent, 1= match, 2=grand */

/******************************************************************
                DEF: Zscore_pq()
set to 1 in clam.h
*****************************************************************/
static int Zscore_pq(ZNODE *node)
{
     if (Zrule_pq==0) return node->exponent;
     else if (Zrule_pq==1) return node->match;
     else return ZZ.matrix[node->i2].grand;
}
/********************************************************************
CAS 10/5/6 Was getting "Fatal error" meaning it could not align bands.
Moreover, it gives difference solutions everytime run, which makes no sense.
Hence, I changed it back to using a constant.
**********************************************************************/
static int get_EXTENT(int c) 
{
int i;
     if (Pz.fp_class==HICF) return HICF_HANG_EXTENT;
     else return AG_HANG_EXTENT;
       // this is not used, since it gives different solutions
       // it performs a wee bit better, so may try to find problem later
     if (Pz.fp_class==HICF) {
          i = GAP(ZZ.matrix[c].nbands);
          if (i < AG_HANG_EXTENT) 
             return  AG_HANG_EXTENT;
          else return i;
     }
     else return AG_HANG_EXTENT;
}
/********************************************************************/
static int get_GAP(int c) 
{
int i;

     if (Pz.fp_class==HICF) return HICF_GAP;
     else return AG_GAP;
      // see note above
     if (Pz.fp_class==HICF) {
          i = GAP(ZZ.matrix[c].nbands);
          if (i < AG_GAP) 
             return  AG_GAP;
          else return i;
     }
     else return AG_GAP;
}

/********************************************************************/
static int get_END_GAP(int w, int c)
{
   if (w==LEFT) return abs(Zset[ZS].leftb - ZZ.matrix[c].leftb);
   else return Zset[ZS].rightb - ZZ.matrix[c].rightb;
}
/*********************************************************************
                    DEF: Zscore_cbmap
called in Zprepare_for_display - which finishs a map
********************************************************************/
static void  Zscore_cbmap()
{
int endb=4;
int i, c, j, g, g1;
int total=0, right=0, left=0;
int p=0, n=0, ends=0, ngap;

    Zset[ZS].n_Qs=0;
    Zset[ZS].score=0.0;
                                         /** CALC SCORES  for CBmap**/
    for (c=0; c< Zset[ZS].n_added_clone; c++) 
    {
        i = Zset[ZS].csort[c].i; 
        ends += MiN(2, ZZ.matrix[i].n_extra);
        total += ZZ.matrix[i].nbands;

        ngap =  get_GAP(ZZ.matrix[i].nbands); 
                    /* count gaps */
        for(g=j=0;j< ZZ.matrix[i].n_cb;j++)
        {
            if (ZZ.matrix[i].cb_index[j+1]!= BIG_NEG_NUM) {
              g1 = abs(ZZ.matrix[i].cb_index[j]-ZZ.matrix[i].cb_index[j+1])-1;
              g += g1;
              if (g1> ngap) g++;
            }
        }
        p+= ZZ.matrix[i].n_cb;

        n+=g;

        g1 = Q_RULE(ZZ.matrix[i].nbands);
        ZZ.matrix[i].cbQ = 0;
             /* no end clones will be Qs */
        if (get_END_GAP(LEFT, i) > ngap && get_END_GAP(RIGHT, i) > ngap) {
            if (g1 <= g+ZZ.matrix[i].n_extra) {
                 if (QTRACE) printf("Set Q: %s gap %d ex %d g1 %d nbands %d \n",
                    arr(acedata, ZZ.matrix[i].cin, CLONE).clone, 
                    g, ZZ.matrix[i].n_extra, g1, ZZ.matrix[i].nbands);
                 ZZ.matrix[i].cbQ = QALIGN;
            } 
        }
        if (ZZ.matrix[i].cbQ == QALIGN) Zset[ZS].n_Qs++;

                    /* find max extras on left and rigth */
        if (ZZ.matrix[i].cbQ!=QALIGN) {
          if (abs(ZZ.matrix[i].leftb - Zset[ZS].leftb) < endb &&
                     abs(ZZ.matrix[i].rightb - Zset[ZS].rightb) < endb) {
             if (left< right) left = MaX(left, ZZ.matrix[i].n_extra);
             else right = MaX(right, ZZ.matrix[i].n_extra);
          }
          else {
            if (abs(ZZ.matrix[i].leftb - Zset[ZS].leftb) < endb)
              left = MaX(left, ZZ.matrix[i].n_extra);
            if (abs(ZZ.matrix[i].rightb - Zset[ZS].rightb) < endb)
              right = MaX(right, ZZ.matrix[i].n_extra);
          }
        }
    }
    left = MaX(0, left-2);
    right = MaX(0, right-2);
    Zset[ZS].score = (float) ((p+left+right) - n) /(float) (total-ends);
    if (Zset[ZS].score < 0) Zset[ZS].score = 0.001;
}
/******************************************************************
                DEF: make_this_Zset()
*****************************************************************/
int make_this_Zset(ZPQZ *pq)
{
ZNODE *node, *Zif_exist();
int i, k, h, b;
int mL=0, hL=0, mR=0, hR=0;
int in_left, in_right;
int fst=-1, last=-1, A=0;
int match, ngap;

    Zrule_avg=tmptest;
    C = pq->C;
    loadCz(&C1z, ZZ.matrix[C].cin);
    Tz.leftb = 1; Tz.rightb = 0;
    match = 0;
    while ((node = Ztraverse(pq->C))!=NULL) {
        if (In_Zset(node->i2))
             match = MaX(match, node->match);
    }

    Tz.high = Q_RULE(match); 
    ngap = get_GAP(C1z.nbands); /* CARI 10.feb.6 ngap replaces all Pz.gap */

    if (Zset[ZS].n_added_clone < 2) {
       for (i=0; i< C1z.nbands; i++) {
            add_extra_ZZ(C, C1z.coords[i], 0);
            Tz.extra[i] = C1z.coords[i];
       }
       Tz.n_extra = i;
       if (Zset[ZS].n_added_clone == 1) {
           maybe_add_Zgroup(C,&mR, &hR, Zset[ZS].hang_right, RIGHT, 0);
           if (mR > 0)
              extend_Zgroup(mR, hR, RIGHT, Zset[ZS].hang_right);
           First_rightb = Zset[ZS].rightb;
           ZZ.matrix[C].rightb = Zset[ZS].rightb;
           ZZ.matrix[C].leftb = Zset[ZS].leftb;
       }
       goto FINISH_ADD;
    }
/* Compare with current bands */
          /* first find window to match */
    if (pq->score == -PQ_BURIED) {
        fst = pq->fstbg = ZZ.matrix[ZZ.matrix[C].high_match].leftb;
        last = pq->lastbg = ZZ.matrix[ZZ.matrix[C].high_match].rightb;
        if (ZZ.matrix[C].mtype==IAM_APPROX) {
          pq->fstbg =  MaX(pq->fstbg - WIN_BURY_FUZZ(C1z.nbands), Zset[ZS].leftb);
          pq->lastbg = MiN(pq->lastbg + WIN_BURY_FUZZ(C1z.nbands), Zset[ZS].rightb);
        }
    }
    else if (Zorder) { 
        for (i=0; i<ZZ.size && ZZ.matrix[i].pick != Zpick-1; i++);
        if (ZZ.matrix[i].pick == Zpick-1) 
           fst = pq->fstbg = ZZ.matrix[i].leftb;
        else {
           fst = pq->fstbg = Zset[ZS].leftb;
           printf("error!!!!!\n");
        }
        last = pq->lastbg = Zset[ZS].rightb;
    }
    else {
       fst = abs(BIG_NEG_NUM);
       last = BIG_NEG_NUM;
       while ((node = Ztraverse(pq->C))!=NULL) {
           if (In_Zset(node->i2)) { 
             fst = MiN(fst, ZZ.matrix[node->i2].leftb);
             last = MaX(last, ZZ.matrix[node->i2].rightb);
           }
       }
       if (last==BIG_NEG_NUM) { /* can happend with Build or Select*/ 
          fst = Zset[ZS].leftb;
          last = Zset[ZS].rightb;
       }
       pq->fstbg =  MaX(fst - WIN_FUZZ_RULE(C1z.nbands), Zset[ZS].leftb);
       pq->lastbg = MiN(last + WIN_FUZZ_RULE(C1z.nbands), Zset[ZS].rightb);
    }
    
              /* get initial extents becaused used in scoring */
    for (i=0; i< C1z.nbands; i++) 
            Tz.extra[i] = C1z.coords[i];
    Tz.n_extra = i;
    if (!Zorder && abs(Zset[ZS].leftb - pq->fstbg) <= ngap)
       maybe_add_Zgroup(C, &mL, &hL, Zset[ZS].hang_left, LEFT, 0);
    if (abs(Zset[ZS].rightb - pq->lastbg) <= ngap)
       maybe_add_Zgroup(C, &mR, &hR,Zset[ZS].hang_right, RIGHT, 0);

      /* only uses mL and mR if a bands are Pz.gap from the end */
      /* it set ZZ.matrix[C].leftb and rightb */
    compare_to_Zset(C, mL, mR, pq->fstbg, pq->lastbg);
   
/* Determine best place */
                 /* Can extend ends ? */
    for (b=i=0; i< C1z.nbands; i++)
        if (Tz.bit[i] == 0) 
            Tz.extra[b++] = C1z.coords[i];
    Tz.n_extra = b;
    mL = mR = 0;

          /* maybe_add_Zgroup only accepts match if > left or right_sp */
    if (!Zorder && (Tz.left_ok ||  !Tz.ok))
        maybe_add_Zgroup(C, &mL, &hL, Zset[ZS].hang_left, LEFT, Tz.left_sp);
    if (Zorder || Tz.right_ok || !Tz.ok) 
        maybe_add_Zgroup(C, &mR, &hR,Zset[ZS].hang_right, RIGHT, Tz.right_sp);

    if (Zadd || Zorder) {
          if (mR>0 && mR>= mL) Tz.right_ok=1;
          else if (mL>0) Tz.left_ok=1;
    }
    if (Tz.right_ok && mR>0 && mR>= mL) h = RIGHT;
    else if (Tz.left_ok && mL>0) h = LEFT;
    else h = 2;

    if (Zset[ZS].n_added_clone <=1) goto GOOD_ADD;

    if (!Zadd && !Zorder && !Tz.ok) {
      if (h == RIGHT)  {
        if (Tz.match+mR < (int) END_RULE(C1z.nbands)) A = ALIGN;
      }
      else if (h == LEFT)  {
          if (Tz.match+mL < (int) END_RULE(C1z.nbands)) A=ALIGN; 
      }
      else A = ALIGN;
    }

    if (A!=ALIGN) goto GOOD_ADD;

    ZZ.matrix[C].cbQ= QALIGN;

    for (i=0; i< Tz.n_extra; i++)
           add_extra_ZZ(C, Tz.extra[i], 0);

    if (Tz.match!=0) { 
      for (i=Tz.leftb, k= 0 ; i<= Tz.rightb; i++, k++) {
        if (Tz.who[k] != -1) 
            add_clone_Zset(i, C, C1z.coords[Tz.who[k]]);
      }
    }
    if (h==LEFT)  {
          maybe_add_Zgroup(C, &mL, &hL, Zset[ZS].hang_left, LEFT, 3);
          if (mL>0) extend_Zgroup(mL, hL, LEFT, Zset[ZS].hang_left);
    }
    else if (h==RIGHT) {
          maybe_add_Zgroup(C, &mR, &hR, Zset[ZS].hang_right, RIGHT, 3);
          if (mR>0) extend_Zgroup(mR, hR, RIGHT, Zset[ZS].hang_right);
    } 
    /* 15dec99 serious heuristics here */
    else {
         maybe_add_Zgroup(C, &mL, &hL, Zset[ZS].hang_left, LEFT, 0);
         maybe_add_Zgroup(C, &mR, &hR, Zset[ZS].hang_right, RIGHT, 0);
         if (mL > mR && mL>8) extend_Zgroup(mL, hL, LEFT, Zset[ZS].hang_left);
         else if (mR>8) extend_Zgroup(mR, hR, RIGHT, Zset[ZS].hang_right);
           /* WN-CAS core dump on Davis short clone stuff */
         else if (Zset[ZS].n_group==0) {
            if (mL > mR && mL>0) extend_Zgroup(mL, hL, LEFT, Zset[ZS].hang_left);
            else if (mR>0) extend_Zgroup(mR, hR, RIGHT, Zset[ZS].hang_right);
            else {fprintf(stderr, "Fatal error in assembly\n");
                 exit(0); }
         }
    }
    msg_clone("2A-", C, Zstep);
    goto FINISH_ADD;

GOOD_ADD: ;
    for (i=0; i< Tz.n_extra; i++)
           add_extra_ZZ(C, Tz.extra[i], 0);

    if (Tz.match!=0) { /* can happen with Zadd */
      in_left = find_Zgroup(Tz.leftb, RIGHT);
      if (in_left==INT_MIN) {
           printf("Bad group in for Tz.leftb \n");
      }        
      in_right = find_Zgroup(Tz.rightb, LEFT);
      if (in_right==INT_MIN) {
           printf("Bad group in for Tz.right \n");
      }        
      for (i=Tz.leftb, k= 0 ; i<= Tz.rightb; i++, k++) {
        if (Tz.who[k] != -1) 
            add_clone_Zset(i, C, C1z.coords[Tz.who[k]]);
        else {
             if ((i > in_left && i < in_right))
                Zset[ZS].band[i].nocnt++; 
        }
      }
    }
    /* do not use Tz.1stg or lastg after roll or extend */
    if (h == LEFT && mL > 0) {
        if (Tz.match!=0) roll_left(Tz.rightb, Tz.lastg, C);
        else ZZ.matrix[C].leftb = Zset[ZS].leftb;
        extend_Zgroup(mL, hL, LEFT, Zset[ZS].hang_left);
    }
    else if (h == RIGHT && mR > 0){ 
        if (Tz.match!=0) roll_right(Tz.leftb, Tz.fstg, C);
        else ZZ.matrix[C].rightb = Zset[ZS].rightb;
        extend_Zgroup(mR, hR, RIGHT, Zset[ZS].hang_right);
    } 
     /* don't roll if group is in the middle, no info on direction to roll */
    else if (Tz.fstg == Tz.lastg) {
        if (Tz.lastg == Zset[ZS].n_group-1) 
             roll_right(Tz.rightb, Tz.lastg, C);
        else if (Tz.fstg == 0)
             roll_left(Tz.leftb, Tz.fstg, C);
    }
    else {
        roll_left(Tz.rightb, Tz.lastg, C);
        roll_right(Tz.leftb, Tz.fstg, C);
    }
    msg_clone("Add", C, Zstep);

FINISH_ADD: ;

    if (Zorder) Zset[ZS].hang_right[0] = C;
    else add_to_hang(C);

    Zset[ZS].n_added_clone++;
    ZZ.matrix[C].pick = Zpick++;

    add_nodes_to_pq(C);

    if (ZZ.matrix[C].leftb > ZZ.matrix[C].rightb) {
        if (Ztrace) printf("**No coords for %s Z(%d, %d) fst(%d,%d)\n",
             arrp(acedata, ZZ.matrix[C].cin, CLONE)->clone, 
             ZZ.matrix[C].leftb, ZZ.matrix[C].rightb, fst,last);
           ZZ.matrix[C].leftb=fst;
           ZZ.matrix[C].rightb=last;
    }
       
    if (Zrule_remove == 1 && Zset[ZS].n_added_clone >=3) 
             remove_any_bad_bands();
return 1;
}

/******************************************************************
                DEF: compare_to_Zset() CMP
*****************************************************************/
int compare_to_Zset(int c, int mL, int mR, int startb, int endb)
{
int  w, i,  k, g, g1, g2, b, b1, b2;
int who,  idiff, save_diff, save_i;
int gt=-1, avg,  match, nomatch=0;
int  Zbit[TBAND];
int loop, goff, lastg,  win_end;
static struct grpy {
   int m, nm, leftb, rightb;
} *grp=NULL;
static int n_grp=0;
int g_save=-1, ok[2], sp[2];
int found, Sok[2], Ssp[2], ngap;
float neg, pos, score, save_score;
void CTG_resources();

   if (startb < Zset[ZS].leftb) startb = Zset[ZS].leftb;
   if (endb > Zset[ZS].rightb) endb = Zset[ZS].rightb;
   
   C = c;
   Tz.leftb = endb+1;
   Tz.rightb = startb-1;
   ngap = get_GAP(C1z.nbands);

   Tz.score= Tz.max=Tz.ok = Tz.match = Tz.nomatch=0;
   Tz.left_ok=Tz.right_ok=Tz.ok=Tz.left_sp=Tz.right_sp=0;
   Tz.right_no = Tz.left_no = 0;
   for (i=0; i<C1z.nbands; i++) Tz.bit[i]=0;
   for (i=0; i<TBAND; i++) Tz.who[i]=-1;

   if (startb > endb) return 0;

   for (i= Zset[ZS].leftb; i <= Zset[ZS].rightb; i++)
         Zset[ZS].band[i].zap = 0;

   if ((endb-startb)+100 > n_win) { /* goes past end if group.right > end */
        n_win = (endb-startb)+200;
        if (win==NULL) win = (struct winy *) malloc((sizeof(struct winy)*n_win));
        else win = (struct winy *) realloc(win, (sizeof(struct winy)*n_win));
        if (win==NULL) {
           CTG_resources();
           exit(1);
        }
        NOMEM3a(win, "win", 0, n_win, sizeof(struct winy)*n_win);
   }

      /* get values for within window to match */
   for (g1=i=0, b=startb; b <= endb && i < n_win-1; )
   {
        g1++;
        g = find_Zgroup(b,2);
        if (g==INT_MIN) {
          printf("Bad group in compare_to_Zset start %d end %d  left %d right %d\n", 
                        startb, endb, Zset[ZS].leftb, Zset[ZS].rightb);
           return 0;
        }        
        while (b <= Zset[ZS].group[g].rightb && i < n_win-1)
        {
           win[i].avg = Zset[ZS].band[b].avg;
           win[i].g = g;
           win[i].who = -1;
           i++; b++;
        }
   }
   win[i].g = -1;
   win_end = i;

   if (g1 > n_grp) { 
        n_grp = g1+100;
        if (grp==NULL) grp =  (struct grpy *) malloc((sizeof(struct grpy)*n_grp));
        else grp = (struct grpy *) realloc(grp, (sizeof(struct grpy)*n_grp));
        NOMEM3a(grp, "grp", 0, n_grp, (sizeof(struct grpy)*n_grp));
   }

   for (i=0; i< TBAND; i++) Zbit[i] = 0;

       /* find matches going from left to right */
   for (match=loop=0; loop <= Zrule_shave; loop++)
     for (w=0; w< win_end && match < C1z.nbands; w++)
     {
           if (win[w].who != -1) continue;
           avg = win[w].avg;
           save_diff =  Zgtol(avg);
           if (loop==1) save_diff *= 2;
           save_i = -1;
           for (k=0; k< C1z.nbands; k++) {
              if (Zbit[k]) continue;
              idiff = abs(avg - C1z.coords[k]);
              if (idiff <= save_diff) {
                  save_diff = idiff;
                  save_i = k;
              }
           }
           if (save_i != -1)
           {
              match++;
              Zbit[save_i] = 1; 
              win[w].who = save_i;
           }
     }
     Tz.max=match;
     if (match==0) goto DONE;

           /* count number of F+ and F- per group */

     for (w=0, goff=g=win[0].g;  w < win_end; g++)
     {
         gt = g-goff;
         if (gt+1 >= n_grp) { /*14may99 Purify array bounds FIX */
            n_grp+=100;
            grp = (struct grpy *) realloc(grp, (sizeof(struct grpy)*n_grp));
            NOMEM3a(grp, "grp", 0,n_grp, sizeof(struct grpy)*n_grp);
         }
         grp[gt].m=grp[gt].nm=0;    
         grp[gt].leftb = w;
         for ( ;win[w].g == g; w++) {
            if (win[w].who== -1) grp[gt].nm++;
            else grp[gt].m++;
         }
         grp[gt].rightb = w-1;
     }
     lastg=gt;
     grp[gt+1].m = 0; 

           /* Score window, slide down one updating grp, win and Zbit structures  
              and score again. Continue until no window. Use best score */
     b1=0;
     while(1)
     {                                
         while (b1 < win_end && win[b1].who == -1) b1++;
         if (b1==win_end) break;

         g1= (win[b1].g - goff);

         for (b2=win_end-1; win[b2].who == -1; b2--);
         if (b1==b2) goto SLIDE;
         g2= (win[b2].g - goff);

         g = g2; save_score=0;
         for (g2=g1; g2 <= g; g2++)
         {
               Ssp[0] = abs(Zset[ZS].leftb - (grp[g1].leftb+startb));
               Ssp[1] = Zset[ZS].rightb - (grp[g2].rightb+startb);
               if (g1!=g2) {
                 Ssp[0] += grp[g1].nm;
                 Ssp[1] += grp[g2].nm;
               }

               Sok[0] = Sok[1] = 0;
               if (Ssp[0] <= ngap && Ssp[1] <= ngap) { 
                   if (mL > mR) Sok[0] = 1;
                   else Sok[1] = 1;
               } 
               else if (Ssp[0] <= ngap) Sok[0]=1;
               else if (Ssp[1] <= ngap) Sok[1]=1;

               if (Sok[0]) score = mL; 
               else if (Sok[1]) score = mR;
               else score = 0;

               if (g1==g2) {
                 pos = grp[g1].m;
                 neg=0;
                 score += pos;
               }
               else {
                 pos = grp[g1].m + grp[g2].m;
                 if (Sok[0]) neg = Ssp[0]; /* grp[g1].nm; */
                 else if (Sok[1]) neg = Ssp[1]; /*grp[g2].nm; */
                 else neg = 0;
                 
                 for (k=g1+1; k<g2; k++) {
                     pos += grp[k].m;
                     neg += (grp[k].nm*1.2);
                 }
                 if (pos > neg/2) score = (score + pos) - neg;
                 else score = 0;
               }
               if (score > save_score || (score == save_score && (int) pos > match)) {
                   save_score = score;
                   g_save = g2;
                   match = pos;
                   nomatch = neg;
                   ok[0] = Sok[0]; ok[1] = Sok[1];
                   sp[0] = Ssp[0]; sp[1] = Ssp[1];
               }
         }
         g2=g_save;
         if(g2>=0)
           for (b2=grp[g2].rightb; b2 > 0 && win[b2].who == -1; b2--);
          
                                
         for (g = g1, k=0; g <= lastg; g++) k += grp[g].m;
         if (k > Tz.max) Tz.max=k;

         if (save_score < Tz.score) goto SLIDE;

         Tz.leftb = startb+b1; Tz.rightb = startb + b2;
         Tz.fstg = win[b1].g; Tz.lastg = win[b2].g;
         Tz.left_sp = sp[0]; Tz.right_sp = sp[1];
         Tz.left_ok = ok[0]; Tz.right_ok = ok[1];
         Tz.match = match;
         Tz.nomatch = nomatch;
         Tz.score = save_score;
         Tz.ok = (match >= Tz.high) ? 1 : 0;

         for (i=0; i< TBAND; i++) { Tz.bit[i] = 0; Tz.who[i] = -1;}
         for (i=0, k=b1; i< TBAND && k<=b2; i++, k++) {
            if (win[k].who != -1) {
               Tz.who[i] = win[k].who;
               Tz.bit[Tz.who[i]]=1;
            }
         }
                             /* slide first match down */
SLIDE: ;
         for (found=loop=0; loop <= Zrule_shave && !found ; loop++)
         {
           for (k=grp[g1].leftb; k<=grp[g1].rightb && !found; k++)
           {
              if (win[k].who== -1) continue;
              who = win[k].who;

              for (i = grp[g1].rightb+1; i < win_end && !found; i++) {
                   if (win[i].who != -1)  continue;

                   idiff =  win[i].avg - C1z.coords[who];
                   if (loop==1) idiff *= 2;
                   if (!Ztol(win[i].avg, idiff)) continue;

                   found=1;
                   win[k].who = -1;
                   grp[g1].m--;
                   grp[g1].nm++;

                   win[i].who = who;
                   g = win[i].g - goff;
                   grp[g].m++;
                   grp[g].nm--;
              }
           } 
         }
         if (!found) {
           for (k=grp[g1].leftb; k<=grp[g1].rightb; k++) {
               if (win[k].who != -1) {
                   Zbit[win[k].who] = 0;
                   win[k].who = -1;
               }
           }
         }
     } 
DONE:;
    if (Zadd || Zorder) {
       if (!Tz.ok && !Tz.left_ok && !Tz.right_ok) {
          if (mR > Tz.match || mL > Tz.match)
             Tz.match = Tz.left_sp = Tz.right_sp = 0;
       }
    }
    if (Tz.match==0) { 
       if (mR>0 && mR > mL) { /* 12dec99 chg check */
          ZZ.matrix[C].leftb = endb;
          ZZ.matrix[C].rightb = endb+1;
          Tz.right_ok = 1;
       }
       else if (mL>0) {
          ZZ.matrix[C].leftb = startb-1;
          ZZ.matrix[C].rightb = startb;
          Tz.left_ok = 1;
       }
       else {
         ZZ.matrix[C].leftb = ((startb+endb)/2) - (ZZ.matrix[C].nbands/2);
         ZZ.matrix[C].rightb = ((startb+endb)/2) + (ZZ.matrix[C].nbands/2);
       }
    }
    else {
       ZZ.matrix[C].rightb = Tz.rightb;
       ZZ.matrix[C].leftb = Tz.leftb;
    }
return Tz.match;
}
/****************************************************************
                  DEF: roll_right
*****************************************************************/
void roll_right(int b, int g, int c)
{
int swap,i;
struct Zset_band tband;

   g = find_Zgroup(b, 2);
   if (g==INT_MIN) {
       printf("Bad group in for roll_right \n");
       return;
   }        
   for (swap=1; swap==1; ) 
   {
       swap = 0;
       for (i=Zset[ZS].group[g].leftb; i<Zset[ZS].group[g].rightb; i++) 
       {
           if (Zset[ZS].band[i].zap==1 && Zset[ZS].band[i+1].zap==0) {
               tband = Zset[ZS].band[i];
               Zset[ZS].band[i] = Zset[ZS].band[i+1];
               Zset[ZS].band[i+1] = tband;
               swap=1;
           }
       }
   }
   i=Zset[ZS].group[g].leftb;  
   if (Zset[ZS].band[i].zap!=1) {
       for (i++; i<=Zset[ZS].group[g].rightb; i++) 
           if (Zset[ZS].band[i].zap==1) {
               if (split_Zgroup(i, g, RIGHT))
                       ZZ.matrix[c].leftb = i;
               break;
            }
   }
}
/****************************************************************
                  DEF: roll_left
*****************************************************************/
void roll_left(int b, int g, int c)
{
int swap,i;
struct Zset_band tband;

   g = find_Zgroup(b, 2);
   if (g==INT_MIN) {
       printf("Bad group in for roll_left \n");
       return;
   }        
   for (swap=1; swap==1; ) 
   {
        swap = 0;
        for (i=Zset[ZS].group[g].rightb; i>Zset[ZS].group[g].leftb; i--) 
        {
            if (Zset[ZS].band[i].zap==1 && Zset[ZS].band[i-1].zap==0) {
                tband = Zset[ZS].band[i];
                Zset[ZS].band[i] = Zset[ZS].band[i-1];
                Zset[ZS].band[i-1] = tband;
                swap=1;
            }
        }
   }
   i=Zset[ZS].group[g].rightb; 
   if (Zset[ZS].band[i].zap!=1) {
       for (i--; i>=Zset[ZS].group[g].leftb; i--) 
          if (Zset[ZS].band[i].zap==1) {
              if (split_Zgroup(i, g, LEFT))
                    ZZ.matrix[c].rightb = i;
              break;
          }
   }
} 
/******************************************************************
                DEF: maybe_add_Zgroup()
*****************************************************************/
void maybe_add_Zgroup(int c, int *M, int *H, int *hang, int w, int limit)
{
int  save_match=0, save_h=-1;
int m,  h,  hh;
int b1, b2, k, l, idiff, lstart, found;

   C = c;
   *M = 0;
   *H = BIG_NEG_NUM;
   save_match = 0;
   if (Tz.n_extra==0) return;
                   /** compare two descending lists **/
   for (h=0; h < HANG; h++)
   {
      hh = hang[h];
      if (hh == -1 || hh == C) continue;
    
      for (lstart = m = k = 0; k < Tz.n_extra; ++k)
      {
          b1 = Tz.extra[k];
          if (b1 == -1) continue;
          for (found=0, l=lstart; l<ZZ.matrix[hh].n_extra && !found; ++l)
          {
              b2 = ZZ.matrix[hh].extra[l];
              if (b2 == -1) continue;

	      idiff = b1 - b2;
              if (Ztol(b1, idiff)) {
	          m++;
	          lstart = l + 1;
	          found=1;
	      }
	      else if (idiff < 0)
	      {
	          lstart = l;
	          found=1;
	      }
           }
       }
       if (m > save_match) {
          if (get_END_GAP(w, hh) < m || Zstep==1) { 
            save_match = m;
            save_h = h;
          }
       }
   }
   if (save_match < limit) return;
   *M = save_match;
   *H = save_h;
}
/****************************************************************
                      DEF: extend_Zgroup
maybe_add must be called before this so C is right
****************************************************************/
void extend_Zgroup(int save_match, int h, int w, int *hang)
{
int G;
int  b, H, nextent;
int b1, b2, k, l, idiff, lstart, found;

   if (save_match == 0) return;
   H = hang[h];
   if ((b = get_Zband(&G, save_match, w)) == BIG_NEG_NUM) return;

   for (lstart = k = 0; k < ZZ.matrix[C].n_extra; ++k)
   {
       b1 = ZZ.matrix[C].extra[k];
       if (b1 == -1) continue;
       for (found=0, l=lstart; l<ZZ.matrix[H].n_extra && !found &&
                               b < Zset[ZS].far_rightb; l++) 
       {
           b2 = ZZ.matrix[H].extra[l];
           if (b2 == -1) continue;

           idiff = b1 - b2;
           if (Ztol(b1, idiff)) {
               add_clone_Zset(b, H, b2);
               add_clone_Zset(b, C, b1);
               ZZ.matrix[C].extra[k] = -1;
               ZZ.matrix[H].extra[l] = -1;

	       lstart = l + 1;
	       found=1;

               b++;
	   }
	   else if (idiff < 0)
	   {
	       lstart = l;
	       found=1;
	   }
        }
    }
    updateZZ(H, w, G);
    updateZZ(C, w, G);

    if (Ztrace>1) printf(" Grp%d (%d, %d) %d-%s  %d-%s match %d\n", G, 
             Zset[ZS].group[G].leftb, Zset[ZS].group[G].rightb,
             C, arr(acedata,ZZ.matrix[C].cin, CLONE).clone, 
             H, arr(acedata,ZZ.matrix[H].cin, CLONE).clone, save_match);

    if (Zadd || Zorder || Zstep==1) return;
                   /** prune extras */
    for (h=0; h < HANG; h++)
    {
        if (hang[h] != -1)  
        {
            nextent = get_EXTENT(ZZ.matrix[hang[h]].nbands); /* CAS 5.10.6 */
            if (w==RIGHT) {
               if (ZZ.matrix[hang[h]].rightb + nextent < Zset[ZS].rightb)
                  hang[h] = -1;
             }
             else {
               if (ZZ.matrix[hang[h]].leftb - nextent > Zset[ZS].leftb)
                  hang[h] = -1;
             }
        }
    }
}
/***********************************************************************
                     DEF: add_to_hang
CARI 10.feb.6  move this from the main routine. 
************************************************************************/
void add_to_hang(int C)
{
int h, nextent;

     nextent = get_EXTENT(ZZ.matrix[C].nbands); 
     if (get_END_GAP(LEFT, C) <= nextent) {
        for (h=0; h<HANG; h++) 
          if (Zset[ZS].hang_left[h]==-1) {
             Zset[ZS].hang_left[h] = C;
             break;
          }
        if (HANG==h) printf("Hit the hange limit\n");
     }
     if (get_END_GAP(RIGHT, C) <= nextent) {
        for (h=0; h<HANG; h++)  
          if (Zset[ZS].hang_right[h]==-1) {
             Zset[ZS].hang_right[h] = C;
             break;
          }
        if (HANG==h) printf("Hit the hange limit\n");
     }
}
/***************************************************************
                  DEF: updataZZ
****************************************************************/
static void updateZZ(int c, int w, int G)
{
int g, e;
int extra[TBAND];

    for (e=g=0; g < ZZ.matrix[c].n_extra; g++)
       if (ZZ.matrix[c].extra[g]!= -1) extra[e++] = ZZ.matrix[c].extra[g];
    for (g=0; g < e; g++) ZZ.matrix[c].extra[g]=extra[g];
    ZZ.matrix[c].extra[g]= -1;
    ZZ.matrix[c].n_extra=e;

    if (w==0) ZZ.matrix[c].leftb = Zset[ZS].group[G].leftb;
    else ZZ.matrix[c].rightb = Zset[ZS].group[G].rightb;
}
/***************************************************************
                  DEF: add_nodes_to_pq
****************************************************************/
void add_nodes_to_pq(int c)
{
ZNODE *node, *Zif_exist();
ZPQZ *pt;
int h;

    C = c;
    while ((node = Ztraverse(C))!=NULL) 
    {
       if (In_Zset(node->i2)) continue;

       if (ZZ.matrix[node->i2].mtype==IAM_EXACT || ZZ.matrix[node->i2].mtype==IAM_APPROX) {
          if (ZZ.matrix[node->i2].high_match  == C) {
             pt = insertZpq(node->i2, PQ_BURIED);
          }
          else { 
             h = ZZ.matrix[node->i2].high_match;
             if (!In_Zset(h)) insertZpq(h, 2);
          }
       }
       else  
       { 
          insertZpq(node->i2,Zscore_pq(node));
       }
    }
}
/***************************************************************
                  DEF: insertZpq
****************************************************************/
ZPQZ *insertZpq(int i1,  int score)
{
int i;

   if (ZZ.matrix[i1].cbmap == ZIGNORED) return NULL;
   if (In_Zset(i1)) return NULL;

   for (i=0; i < n_Zpq; i++) 
      if (Zpq[i].C == i1 && Zpq[i].score != -1) {
          if (Zrule_pq<2) Zpq[i].score = MaX(Zpq[i].score, score); 
          return &(Zpq[i]);
      }
   for (i=0; i < n_Zpq; i++) 
      if (Zpq[i].score < 0) break;
      
   if (i == n_Zpq) { 
      if (n_Zpq == max_Zpq) {
         max_Zpq += 30;
         if (Zpq==NULL) Zpq = (ZPQZ *) calloc(max_Zpq, sizeof(ZPQZ));
         else Zpq = (ZPQZ *) realloc(Zpq, max_Zpq * sizeof(ZPQZ));
         NOMEM3a(Zpq,"Zpq", NULL, max_Zpq, max_Zpq * sizeof(ZPQZ));
         for (i=n_Zpq; i<max_Zpq; i++) {
            Zpq[i].C = Zpq[i].fstbg = Zpq[i].lastbg = Zpq[i].score = 0;
            Zpq[i].score = -1;
         }
      }
      i = n_Zpq;
      n_Zpq++;
   }
   Zpq[i].C = i1;
   Zpq[i].score = score;
   return &(Zpq[i]);
}
/***************************************************************
                  DEF: popZpq
****************************************************************/
ZPQZ *popZpq()
{
ZNODE *node, *Zif_exist();
int i, j, score;

POP: ;
   if (n_Zpq==0) return NULL;
   for (j=-1, score= 0, i=0; i < n_Zpq; i++) 
   {
      if (Zpq[i].score > score ) {
         score = Zpq[i].score;
         j = i;
      }
   }
   if (j== -1) return NULL;

      /* parent can pop before any peers have been aligned */
   if (Zpq[j].score==2) {
      while ((node = Ztraverse(Zpq[j].C))!=NULL) 
      {
        if (ZZ.matrix[node->i2].mtype==IAM_EXACT || ZZ.matrix[node->i2].mtype==IAM_APPROX) 
        {
          if (ZZ.matrix[node->i2].high_match  == Zpq[j].C) {
             ZZ.matrix[node->i2].high_match  = -1;
             ZZ.matrix[node->i2].mtype  = IAM_NOTHING;
             insertZpq(node->i2,Zscore_pq(node));
          }
        }
      }
      ZZ.matrix[Zpq[j].C].mtype  = IAM_NOTHING;
      Zpq[j].score = 3;
      goto POP; 
   }
   if (Zpq[j].score !=  PQ_BURIED) Zpq[j].score = -1;
   else Zpq[j].score = -PQ_BURIED;
   return &(Zpq[j]);
}
/***************************************************************
                  DEF: drainZpq
****************************************************************/
ZPQZ *drainZpq()
{
int i, j, score;

   if (n_Zpq==0) return NULL;
   for (j=-1, score= 0, i=0; i < n_Zpq; i++) 
   {
      if (Zpq[i].score > -1 ) {
         j = i;
         break;
      }
   }
   if (j== -1) return NULL;
   Zpq[j].score = -1;
   return &(Zpq[j]);
}
/*******************************************************
                    DEF: add_clone_Zset
********************************************************/
void add_clone_Zset(int  b, int c, int cband)
{
int i, j, k, cnt;

   if (b < Zset[ZS].leftb || b > Zset[ZS].rightb) {
      BAD("FPC error: add_clone_Zset");
      exit(0);
   }
   if (Zset[ZS].band[b].cnt == Zset[ZS].band[b].max) {
      Zset[ZS].band[b].max+= 20;
      i = Zset[ZS].band[b].max * sizeof(int);
      if (Zset[ZS].band[b].clone == NULL) {
         Zset[ZS].band[b].clone = (int *) malloc(i);
         Zset[ZS].band[b].val = (int *) malloc(i);
      }
      else {
         Zset[ZS].band[b].clone = (int *) realloc(Zset[ZS].band[b].clone, i);
         Zset[ZS].band[b].val = (int *) realloc(Zset[ZS].band[b].val, i);
      }
      NOMEM2a(Zset[ZS].band[b].clone,"band-clone", Zset[ZS].band[b].max, i);
      NOMEM2a(Zset[ZS].band[b].val,"band-vaL",Zset[ZS].band[b].max, i);
   }

   i = Zset[ZS].band[b].cnt;
   Zset[ZS].band[b].cnt++;
   Zset[ZS].band[b].clone[i] = c;
   Zset[ZS].band[b].val[i] = cband;
   Zset[ZS].band[b].zap = 1;
   Zset[ZS].band[b].total += cband;
   Zset[ZS].band[b].high = MaX(cband, Zset[ZS].band[b].high);
   Zset[ZS].band[b].low = MiN(cband, Zset[ZS].band[b].low);
    
   if (Zrule_avg==0) { /* average */
      if (Zset[ZS].band[b].cnt!=0) 
         Zset[ZS].band[b].avg = Zset[ZS].band[b].total / Zset[ZS].band[b].cnt;
      else {
         fprintf(stderr,"Internal clam error: Band cnt = 0\n");
         Zset[ZS].band[b].avg = 1;
      }
   }
   else if (Zrule_avg==1) { /* median */
      cnt = Zset[ZS].band[b].cnt; 
      for (i=0; i<cnt-1; i++)
         for (j=i+1; j<cnt; j++)
            if (Zset[ZS].band[b].val[i] > Zset[ZS].band[b].val[j]) 
            {
                k = Zset[ZS].band[b].val[i]; 
                Zset[ZS].band[b].val[i] = Zset[ZS].band[b].val[j]; 
                Zset[ZS].band[b].val[j] = k; 
                k = Zset[ZS].band[b].clone[i]; 
                Zset[ZS].band[b].clone[i] = Zset[ZS].band[b].clone[j]; 
                Zset[ZS].band[b].clone[j] = k; 
            }
      i = cnt >> 1;
      if (cnt % 2 == 0) 
        Zset[ZS].band[b].avg = 
         (Zset[ZS].band[b].val[i-1] +  Zset[ZS].band[b].val[i]) >> 1;
      else 
        Zset[ZS].band[b].avg = Zset[ZS].band[b].val[i];
   }
   else if (Zrule_avg==2) /* midpt */
      Zset[ZS].band[b].avg = (Zset[ZS].band[b].high +  Zset[ZS].band[b].low) >> 1;
 
}
/*******************************************************
                    DEF: add_extra_ZZ
********************************************************/
void add_extra_ZZ(int c, int  b, int sort)
{
int i, j, k;

   if (ZZ.matrix[c].n_extra == ZZ.matrix[c].max_extra) {
      ZZ.matrix[c].max_extra += 20;
      i = ZZ.matrix[c].max_extra*sizeof(int);
      if (ZZ.matrix[c].extra==NULL) ZZ.matrix[c].extra = (int *) malloc(i);
      else ZZ.matrix[c].extra = (int *) realloc(ZZ.matrix[c].extra, i);
      NOMEM2a(ZZ.matrix[c].extra,"extra", ZZ.matrix[c].max_extra,i);
   }
   i= ZZ.matrix[c].n_extra;
   ZZ.matrix[c].n_extra++;
   ZZ.matrix[c].extra[i] = b;
   if (!sort) return;

   for (i= 0; i<ZZ.matrix[c].n_extra-1; i++)
      for (j= i+1; j<ZZ.matrix[c].n_extra; j++) 
         if (ZZ.matrix[c].extra[i] > ZZ.matrix[c].extra[j]) {
             k = ZZ.matrix[c].extra[i];
             ZZ.matrix[c].extra[i] = ZZ.matrix[c].extra[j];
             ZZ.matrix[c].extra[j] = k;
         }
}
/********************************************************
                    DEF: find_Zgroup
could be binary seach
********************************************************/
int find_Zgroup(int b, int w)
{ 
int i;

  for (i=0; i< Zset[ZS].n_group; i++)
  {
    if (b >= Zset[ZS].group[i].leftb &&
        b <= Zset[ZS].group[i].rightb) 
    {
             if (w==LEFT) return Zset[ZS].group[i].leftb;
             else if (w==RIGHT) return Zset[ZS].group[i].rightb;
             else return i;
    }
  }
  printf("CB error: cannot find group for %d %s %d %d %d - continuing.. \n",
         b, C1z.clone, Zstep, ZS, Zset[ZS].n_added_clone);
return INT_MIN;
}
/********************************************************
                    DEF: split_Zgroup
index where element should be;
********************************************************/
static int split_Zgroup(int i, int g, int w)
{ 
int j;

  if (g==INT_MIN) return 0;
  if (Zset[ZS].group[g].leftb == Zset[ZS].group[g].rightb) return 0;
  if (w==LEFT && i+1 == Zset[ZS].group[g].leftb) return 0;
  if (w==RIGHT && i-1 == Zset[ZS].group[g].rightb) return 0;

  if (Ztrace>1) printf("Split%d %d %d, %d at %d\n", w, 
      g, Zset[ZS].group[g].leftb, Zset[ZS].group[g].rightb, i);

  if (w==LEFT) {
      if (Zrule_no_1_grp) {
         if (abs(Zset[ZS].group[g].leftb - i) == 1) return 0;
         if (abs(Zset[ZS].group[g].rightb - (i+1)) == 1) return 0;
      }
      if (i+1 == Zset[ZS].group[g].leftb) return 0;
      j = Zset[ZS].group[g].rightb;
      Zset[ZS].group[g].rightb = i;
      get_Zgroup(i+1, j);
  }
  else {
      if (Zrule_no_1_grp) {
          if (abs(Zset[ZS].group[g].rightb - i) == 1) return 0;
          if (abs(Zset[ZS].group[g].leftb - (i-1)) == 1) return 0;
      }
      if (i-1 == Zset[ZS].group[g].rightb) return 0;
      j = Zset[ZS].group[g].leftb;
      Zset[ZS].group[g].leftb = i;
      get_Zgroup(j, i-1);
  }
return 1;
}
/*******************************************************
                    DEF: get_Zgroup
********************************************************/
int get_Zgroup(int l, int r)
{
int g, i, j;
struct Zset_group tmp;

   if (Zset[ZS].n_group == Zset[ZS].max_group) {
       Zset[ZS].max_group += 50;
       i=Zset[ZS].max_group*sizeof(struct Zset_group);
       if (Zset[ZS].group==NULL) Zset[ZS].group = (struct Zset_group *) malloc(i);
       else Zset[ZS].group = (struct Zset_group *) realloc(Zset[ZS].group, i);
       NOMEM3a(Zset[ZS].group, "Zgroup", -1, Zset[ZS].max_group,i);
       for (i=Zset[ZS].n_group; i < Zset[ZS].max_group; i++) {
          Zset[ZS].group[i].leftb = 0;
          Zset[ZS].group[i].rightb = 0;
       } 
   }
   g = Zset[ZS].n_group;
   Zset[ZS].n_group++;
   Zset[ZS].group[g].leftb = l;
   Zset[ZS].group[g].rightb = r;

   for (i=0; i< Zset[ZS].n_group-1; i++)
     for (j=i+1; j< Zset[ZS].n_group; j++)
        if (Zset[ZS].group[i].leftb > Zset[ZS].group[j].leftb) {
           tmp = Zset[ZS].group[j];
           Zset[ZS].group[j] = Zset[ZS].group[i];
           Zset[ZS].group[i] = tmp;
        }
return 1;
}
/*******************************************************
                    DEF: init_bands
********************************************************/
static void init_bands(int a, int b) {
int i;
     for (i=a; i < b; i++) {
         Zset[ZS].root[i].high = BIG_NEG_NUM;
         Zset[ZS].root[i].low = 99999;
         Zset[ZS].root[i].avg = Zset[ZS].root[i].total = 0;
         Zset[ZS].root[i].max = Zset[ZS].root[i].cnt = Zset[ZS].root[i].nocnt = 0;
         Zset[ZS].root[i].clone = NULL;
         Zset[ZS].root[i].val = NULL;
     } 
}
/*******************************************************
                    DEF: get_Zband
********************************************************/
static int get_Zband(int *G, int m, int w)
{
int b, i, k;
int inc=10, add, old_root_right;

   if (m==0) {
     printf("CB error: alloc 0 bands\n");
     Zprepare_for_display(); Zclamdump();
     return 0;
   }
   if (w == 0)
   {
     if (abs(Zset[ZS].far_leftb) - abs(Zset[ZS].leftb) <= m) {
         add = m+inc;
         Zset[ZS].max_band += add;
         i = sizeof(struct Zset_band)*Zset[ZS].max_band;
         if (Zset[ZS].root==NULL) Zset[ZS].root = (struct Zset_band *) malloc(i);
         else Zset[ZS].root = (struct Zset_band *) realloc(Zset[ZS].root, i);
         NOMEM3a(Zset[ZS].root, "Zroot-L", BIG_NEG_NUM,Zset[ZS].max_band,i);

         old_root_right =  abs(Zset[ZS].far_leftb) + Zset[ZS].rightb;
         for (i= old_root_right+add, k = old_root_right; k>=0; i--, k--)
            Zset[ZS].root[i] = Zset[ZS].root[k];
         init_bands(0, add);
         init_bands(old_root_right+add+1, Zset[ZS].max_band);
         Zset[ZS].far_leftb -= (m+inc);
      }
      if (Zrule_no_1_grp && m==1) Zset[ZS].group[0].leftb--;
      else get_Zgroup(Zset[ZS].leftb - m, Zset[ZS].leftb-1);
      *G = 0;
      Zset[ZS].leftb -=  m;
      b = Zset[ZS].leftb;
   }
   else {
      if (Zset[ZS].far_rightb - Zset[ZS].rightb <= m) {
         Zset[ZS].max_band += m+inc;
         i = sizeof(struct Zset_band)*Zset[ZS].max_band;
         if (Zset[ZS].root==NULL) Zset[ZS].root = (struct Zset_band *) malloc(i);
         else Zset[ZS].root = (struct Zset_band *) realloc(Zset[ZS].root, i);
         NOMEM3a(Zset[ZS].root, "Zroot-R", BIG_NEG_NUM,Zset[ZS].max_band,i);
         init_bands(abs(Zset[ZS].far_leftb)+Zset[ZS].rightb+1,Zset[ZS].max_band);
         Zset[ZS].far_rightb += m+inc;
      } 
      if (Zrule_no_1_grp && m==1) Zset[ZS].group[Zset[ZS].n_group-1].rightb++;
      else  get_Zgroup(Zset[ZS].rightb+1, Zset[ZS].rightb+m);
      *G = Zset[ZS].n_group-1;
      b = Zset[ZS].rightb+1;
      Zset[ZS].rightb += m;
   }
   Zset[ZS].band = Zset[ZS].root - Zset[ZS].far_leftb;
   Zset[ZS].n_band += m;
return b;
}
/*****************************************************************/
static void msg_clone(char *msg, int C, int z)
{
if (Ztrace==0) return;
   printf(
" %s%d %d-%s  score %d match %d (high %d) no %d  left %d/%d right %d/%d ex %d max %d %d\n",
        msg, z, C, arr(acedata,ZZ.matrix[C].cin, CLONE).clone, 
        Tz.score, Tz.match, Tz.high, Tz.nomatch, Tz.leftb, Zset[ZS].leftb,  
        Tz.rightb, Zset[ZS].rightb, Tz.n_extra, Tz.max, ZZ.matrix[C].high_match);
}
/*************************************************************
                       DEF: remove_any_bad_bands()
**************************************************************/
int remove_any_bad_bands()
{
int i, j, b, c, t, g, e;
struct Zset_band *tmp;
struct Zset_group gtmp;
int del=0, found=0, first,  right;

    for (b=Zset[ZS].leftb; b <= Zset[ZS].rightb && !found; b++)
       if (Zset[ZS].band[b].cnt <= Zset[ZS].band[b].nocnt && 
               Zset[ZS].band[b].nocnt> 2) found=1; 
       else if (Zset[ZS].band[b].cnt <=  1) found=1;
    if (!found) return 0;

                       /** find this and rest of bad bands **/

    tmp = (struct Zset_band *) calloc(Zset[ZS].n_band,sizeof(struct Zset_band));
    NOMEM3(tmp, "remove bands", 0);
    first=b-1;
    right = Zset[ZS].rightb;

    for (t=0, b--; b <= right; b++)
    {
       found=0;
       if (Zset[ZS].band[b].cnt <= Zset[ZS].band[b].nocnt && 
               Zset[ZS].band[b].nocnt> 2) found=1; 
       else if (Zset[ZS].band[b].cnt <=  1) found=1;
       if (!found)
       { 
          tmp[t++] = Zset[ZS].band[b];
          continue;
       } 
       if (Ztrace>1) printf("remove %d (%d, %d)\n",b,
            Zset[ZS].band[b].cnt, Zset[ZS].band[b].nocnt);
       remove_bad++;
                                            /** find band and add to extra **/
       for (i=0; i < Zset[ZS].band[b].cnt; i++)
       {   c = Zset[ZS].band[b].clone[i];
           e = Zset[ZS].band[b].val[i];
           loadCz(&C1z, ZZ.matrix[c].cin);
           for (j=0; j< C1z.nbands; j++)
               if (C1z.coords[j] == e) {
                    add_extra_ZZ(c, C1z.coords[j], 1);
                    break;
               }
           if (j == C1z.nbands) {
               printf("XX %s e %d\n", C1z.clone, e);
               printf("ZS %d band %d index %d clone %d val %d cnt %d\n",
                    ZS, b, i, c, e, Zset[ZS].band[b].cnt);
           }
        }

                                           /** remove from group **/
        g = find_Zgroup(b-del, 2);
        if (g==INT_MIN) {
             printf("Bad group in  remove_band_bands\n");
             continue;
        }        
        for (i=g+1; i < Zset[ZS].n_group; i++) {
           Zset[ZS].group[i].leftb--;
           Zset[ZS].group[i].rightb--;
        }
        Zset[ZS].group[g].rightb--;
        if (Zset[ZS].group[g].rightb < Zset[ZS].group[g].leftb) 
        {
              Zset[ZS].group[g].leftb = 99999;
              for (i=0; i< Zset[ZS].n_group-1; i++)
                for (j=i+1; j< Zset[ZS].n_group; j++)
                   if (Zset[ZS].group[i].leftb > Zset[ZS].group[j].leftb) {
                      gtmp = Zset[ZS].group[j];
                      Zset[ZS].group[j] = Zset[ZS].group[i];
                      Zset[ZS].group[i] = gtmp;
                   }
              Zset[ZS].n_group--; 
        }
                                           /** adjust clone coords **/
        for (i=0; i< ZZ.size; i++)
           if (ZZ.matrix[i].cbmap == ZS && In_Zset(i)) { /* cari 17feb05 */
               if (ZZ.matrix[i].rightb >= b) ZZ.matrix[i].rightb--;
               if (ZZ.matrix[i].leftb > b) ZZ.matrix[i].leftb--;
           }
        Zset[ZS].n_band--;
        Zset[ZS].rightb--;
        del++;
        ZFREE(Zset[ZS].band[b].clone,"band[b].clone"); 
        ZFREE(Zset[ZS].band[b].val,"band[b].val"); 
    }
    for (i=first, j=0; j < t; i++, j++) Zset[ZS].band[i] = tmp[j];
    for (; i <= right; i++) {
         Zset[ZS].band[i].high = BIG_NEG_NUM;
         Zset[ZS].band[i].low = 99999;
         Zset[ZS].band[i].avg = Zset[ZS].band[i].total = 0;
         Zset[ZS].band[i].max = Zset[ZS].band[i].cnt = Zset[ZS].band[i].nocnt = 0;
         Zset[ZS].band[i].clone = NULL;
         Zset[ZS].band[i].val = NULL;
     } 
    free(tmp);
                        /** delete any group if necessay **/
    for (i=0; i< Zset[ZS].n_group-1; i++)
      for (j=i+1; j< Zset[ZS].n_group; j++)
        if (Zset[ZS].group[i].leftb > Zset[ZS].group[j].leftb) {
           gtmp = Zset[ZS].group[j];
           Zset[ZS].group[j] = Zset[ZS].group[i];
           Zset[ZS].group[i] = gtmp;
        }
return del;
}
/********************************************************************
                   DEF: add_extra_band
***********************************************************************/
static void add_extra_band(int k, int c, int b)
{
int x;

      add_clone_Zset(k, c, b);
      x = ZZ.matrix[c].n_cb;
      ZZ.matrix[c].cb_index[x] = k;
      ZZ.matrix[c].cb_index[x+1] = BIG_NEG_NUM;
      ZZ.matrix[c].n_cb++;
}

/************************************************************
                DEF: Zprepare_for_display
*************************************************************/
void Zprepare_for_display()
{
int i, j, c, k, x, e;
int left=0, right=0, chg, extra[TBAND], last, cnt=0, cnt2=0, big, found;
int n_cb=0;
int save_Zsort;

/* make clone who list for draw routine */
    if (Zset[ZS].n_group == 0) return;

    save_Zsort = Zsort;
    Zsort = 0;
    cbmap_sort(); /* cb_display.c - sorts by order they were picked */
    Zsort = save_Zsort;

    if (Zadd)
      for (i=0; i<ZZ.size; i++) ZZ.matrix[i].n_cb=0;

                                         /** assign bands to clones **/
    for (i=Zset[ZS].leftb; i<= Zset[ZS].rightb; i++)
    { 
       for (j=0; j < Zset[ZS].band[i].cnt; j++) 
       {
           c = Zset[ZS].band[i].clone[j];
           if (ZZ.matrix[c].cbmap != ZS || ZZ.matrix[c].cb_index==NULL) {
               printf("CB error: inconsistent band to clone %d %d\n",i,c);
               printf("band cnt %d\n",Zset[ZS].band[i].cnt);
           }
           else {
             k=ZZ.matrix[c].n_cb;
             ZZ.matrix[c].cb_index[k] = i;
             ZZ.matrix[c].cb_index[k+1] = BIG_NEG_NUM;
             ZZ.matrix[c].n_cb++;
           } 
       }
    }
             /* the average moves around so much that quite a few can be added now,
                even though they were not in the range when tested earlier */

    for (i=0; i<Zset[ZS].n_added_clone; i++) 
    {
         c = Zset[ZS].csort[i].i;
         for (j=0; j<ZZ.matrix[c].n_extra; j++)
           extra[j] = ZZ.matrix[c].extra[j];
         
         last = ZZ.matrix[c].cb_index[0];
         n_cb = ZZ.matrix[c].n_cb;
         for(chg=0, j=1; j < n_cb;j++)
         {
                      /* if two bands are not consecutive, can try to fill in */
              for (k=last+1; k<ZZ.matrix[c].cb_index[j]; k++)
              {
                   if (k< BIG_NEG_NUM+1) {
                      printf("%d %d %d %d\n",k, 
                           ZZ.matrix[c].cb_index[j-1], n_cb, ZZ.matrix[c].n_cb);
                      break;
                   }
                   band = Zset[ZS].band[k].avg;
                   for (found=e=0; !found && e<ZZ.matrix[c].n_extra ; e++)
                   {
                        if (extra[e]== -1) continue;
                        big = Zgtol(band);
                        if (abs(band - extra[e]) <= big << 1) {
                               cnt++;
                               chg=found=1;
                               add_extra_band(k, c, extra[e]);
                               extra[e] = -1;
                        }
    
                    }
               }
               last = ZZ.matrix[c].cb_index[j];
         }
                                   /* exterior */ 
         left=0;
         last=1; 
         j=ZZ.matrix[c].leftb-1; 
         k=ZZ.matrix[c].rightb+1;

         while (left==0 || right==0)
         {
              if (last==0 && !left){
                  x = j;
                  j--;
                  last = 1; 
                  if (x < Zset[ZS].leftb || abs(x - ZZ.matrix[c].leftb) > 2) {
                      left=1;
                      continue;
                  }
              }
              else {
                  x = k;
                  k++;
                  last = 0;
                  if (x > Zset[ZS].rightb ||  abs(x - ZZ.matrix[c].rightb) > 2) {
                      right=1;
                      continue;
                  }
              }
              band = Zset[ZS].band[x].avg;
              for (e=found=0; !found && e<ZZ.matrix[c].n_extra ; e++)
              {
                    if (extra[e]== -1) continue;
                    big = Zgtol(band);
                    if (abs(band - extra[e]) <= big << 1) {
                        chg = found = 1;
                        cnt2++;
                        add_extra_band(x, c, extra[e]);
                        if (last==1) ZZ.matrix[c].leftb = x;
                        else ZZ.matrix[c].rightb = x;
                        extra[e] = -1;
                    }
              }
         }

         if (!chg) continue;

         for (k=0; k< ZZ.matrix[c].n_cb; k++)
           for (j=k+1; j< ZZ.matrix[c].n_cb; j++)
               if (ZZ.matrix[c].cb_index[k] > ZZ.matrix[c].cb_index[j]) {
                 x = ZZ.matrix[c].cb_index[k];
                 ZZ.matrix[c].cb_index[k] = ZZ.matrix[c].cb_index[j];
                 ZZ.matrix[c].cb_index[j] = x;
               }
          ZZ.matrix[c].cb_index[k] = BIG_NEG_NUM;

          x = ZZ.matrix[c].n_extra;
          ZZ.matrix[c].n_extra=0;
          for (k=e=0; e< x; e++) {
              if (extra[e]!= -1) {
                 ZZ.matrix[c].n_extra++;
                 ZZ.matrix[c].extra[k++] = extra[e];
              }
          }
    }
    Zscore_cbmap();
}
